Historical Renewables Era#
This tutorial explores the Historical_2023_etys scenario using the full ETYS network topology (2000+ buses). By 2023, Great Britain’s electricity system had undergone significant transformation with coal phase-out complete and wind/solar providing substantial generation.
What You’ll Learn#
Working with the detailed ETYS transmission network
Analyzing high renewable penetration scenarios
Identifying transmission constraints and congestion hotspots
Understanding renewable curtailment patterns
Evaluating system flexibility requirements
Prerequisites#
Run the workflow to generate the solved network:
snakemake resources/network/Historical_2023_etys_solved.nc -j 4
Scenario Overview#
Parameter |
Value |
|---|---|
Modelled Year |
2023 |
Network Model |
ETYS (2000+ buses) |
Data Source |
DUKES (thermal), REPD (renewables), ESPENI (demand) |
Solve Period |
First week of January |
Key Features |
No coal, high wind/solar, network detail |
[1]:
# Import required libraries
import os
import sys
import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
from plotly.subplots import make_subplots
from pyproj import Transformer
from pathlib import Path
# Configure matplotlib and plotly
plt.rcParams.update({'font.size': 12})
plt.style.use('ggplot')
%matplotlib inline
# Configure Plotly to output HTML for Sphinx/nbsphinx compatibility
pio.renderers.default = 'notebook_connected'
# Suppress warnings
import warnings
warnings.filterwarnings('ignore')
print('✓ Libraries imported successfully')
✓ Libraries imported successfully
[2]:
# Load solved network (resolve path robustly from different working directories)
network_rel = Path('resources/network/Historical_2023_etys_solved.nc')
network_path = network_rel
# If the relative path doesn't exist from the current cwd, search upward for the repository root
if not network_path.exists():
for parent in [Path.cwd()] + list(Path.cwd().parents)[:6]:
candidate = parent / network_rel
if candidate.exists():
network_path = candidate
break
else:
# leave as relative; we'll raise a clear error below if not found
network_path = network_rel
print(f'Loading network from: {network_path}')
if not network_path.exists():
raise FileNotFoundError(f'Network file not found. Tried: {network_path}.\\nMake sure "resources/network/Historical_2023_etys_solved.nc" exists relative to the repository root.')
network = pypsa.Network(str(network_path))
# Derive a scenario id from the network filename when not provided
scenario_id = 'Historical_2023_etys'
print(f'✓ Network loaded successfully')
print(f' Buses: {len(network.buses)}')
print(f' Generators: {len(network.generators)}')
print(f' Lines: {len(network.lines)}')
print(f' Storage units: {len(network.storage_units)}')
print(f' Timesteps: {len(network.snapshots)}')
print(f' Period: {network.snapshots[0]} to {network.snapshots[-1]}')
Loading network from: c:\Users\alyden\OneDrive - University of Edinburgh\Python\PyPSA-GB v0.0.1\resources\network\Historical_2023_etys_solved.nc
INFO:pypsa.network.io:Imported network 'Historical_2023_etys (Full)' has buses, carriers, generators, lines, links, loads, storage_units, sub_networks, transformers
✓ Network loaded successfully
Buses: 2044
Generators: 4766
Lines: 1592
Storage units: 81
Timesteps: 168
Period: 2023-05-15 00:00:00 to 2023-05-21 23:00:00
Network Summary#
[3]:
# Calculate generation by carrier
p_by_carrier = network.generators_t.p.groupby(
network.generators.carrier, axis=1).sum()
# Add storage output
storage_by_carrier = network.storage_units_t.p.groupby(
network.storage_units.carrier, axis=1).sum()
storage_by_carrier[storage_by_carrier < 0] = 0
p_by_carrier = pd.concat([p_by_carrier, storage_by_carrier], axis=1)
# Add H2_turbine links (hydrogen power plants modeled as links, not generators)
if len(network.links) > 0:
h2_turbine_links = network.links[network.links['carrier'] == 'H2_turbine']
if len(h2_turbine_links) > 0 and len(network.links_t.p1) > 0:
# p1 is the power output at bus1 (GB side), which is negative
h2_output = network.links_t.p1[h2_turbine_links.index].abs()
h2_ts = pd.DataFrame({'H2_turbine': h2_output.sum(axis=1)})
p_by_carrier = pd.concat([p_by_carrier, h2_ts], axis=1)
print(f'Added H2_turbine generation: {h2_turbine_links.p_nom.sum():.0f} MW capacity')
# Separate interconnectors from internal power lines
# Interconnectors connect to external buses (e.g., 'HVDC_External_France')
if len(network.links) > 0:
# Identify interconnectors vs internal links by checking for 'external' in bus names
is_interconnector = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
interconnector_links = network.links[is_interconnector].index
internal_links = network.links[~is_interconnector].index
# Get interconnector imports (positive p0 = flow into GB from external)
if len(interconnector_links) > 0:
ic_flows = network.links_t.p0[interconnector_links].copy()
ic_flows[ic_flows < 0] = 0
interconnector_import = pd.DataFrame({'Interconnectors Import': ic_flows.sum(axis=1)})
p_by_carrier = pd.concat([p_by_carrier, interconnector_import], axis=1)
# Report internal links separately (not included in generation mix)
if len(internal_links) > 0:
print(f'Note: {len(internal_links)} internal power transmission links (not shown in generation mix)')
else:
print('No links in network')
print('Generation Mix (MWh):')
print(p_by_carrier.sum().sort_values(ascending=False))
print()
print(f'Total demand satisfied: {network.loads_t.p_set.sum().sum():,.0f} MWh')
Note: 3 internal power transmission links (not shown in generation mix)
Generation Mix (MWh):
wind_offshore 1.150376e+06
EU_import 9.750260e+05
wind_onshore 6.977331e+05
nuclear 5.686804e+05
CCGT 4.091246e+05
waste_to_energy 3.105649e+05
solar_pv 2.760712e+05
coal 1.719459e+05
Interconnectors Import 1.631320e+05
landfill_gas 1.024287e+05
Pumped Storage Hydroelectricity 4.560986e+04
biogas 4.497526e+04
advanced_biofuel 4.392902e+04
OCGT 3.257013e+04
large_hydro 2.311478e+04
Battery 2.195361e+04
sewage_gas 7.551726e+03
small_hydro 5.620384e+03
shoreline_wave 2.765148e+03
tidal_stream 6.709302e+01
load_shedding 6.912225e+00
dtype: float64
Total demand satisfied: 4,805,482 MWh
Interactive Generation Mix Visualization#
[4]:
# Interactive stacked area chart with Plotly
# Define distinct colors for different technologies (ordered: baseload -> renewables -> peakers)
color_map = {
# Baseload (bottom of stack)
'nuclear': '#9467BD', # Purple - distinctive for nuclear
# Low-carbon dispatchable
'waste_to_energy': '#8C564B', # Brown
'landfill_gas': '#D2691E', # Chocolate
'biogas': '#2E8B57', # Sea green
'sewage_gas': '#556B2F', # Dark olive
'advanced_biofuel': '#6B8E23', # Olive drab
'biomass': '#228B22', # Forest green
'Other Bioenergy': '#3CB371', # Medium sea green (grouped)
# Renewables (variable)
'wind_offshore': '#1E90FF', # Dodger blue
'wind_onshore': '#00CED1', # Dark turquoise
'solar_pv': '#FFD700', # Gold
'large_hydro': '#4169E1', # Royal blue
'small_hydro': '#87CEEB', # Sky blue
'tidal_stream': '#20B2AA', # Light sea green
'shoreline_wave': '#48D1CC', # Medium turquoise
'Other Renewables': '#5F9EA0', # Cadet blue (grouped)
# Storage (mid-merit)
'battery': '#32CD32', # Lime green
'Pumped Storage Hydroelectricity': '#00FA9A', # Medium spring green
'pumped_hydro': '#00FA9A',
# Gas (flexible)
'CCGT': '#FF6347', # Tomato red
'OCGT': '#FF4500', # Orange red (peaker)
# Coal and other fossil
'coal': '#696969', # Dim grey
# Imports
'Interconnectors Import': '#DDA0DD', # Plum
'EU_import': '#DA70D6', # Orchid
# Hydrogen system
'H2': '#00FFFF', # Cyan (legacy H2 generators)
'H2_turbine': '#00CED1', # Dark turquoise (H2 to power)
'electrolysis': '#40E0D0', # Turquoise (power to H2)
'H2_gas': '#7FFFD4', # Aquamarine (H2 storage)
# Emergency
'load_shedding': '#DC143C' # Crimson
}
# Define stacking order: baseload at bottom, peakers at top
stacking_order = [
'nuclear', # Baseload (bottom)
'waste_to_energy', 'landfill_gas', 'biogas', 'sewage_gas', 'advanced_biofuel', 'biomass', 'Other Bioenergy',
'wind_offshore', 'wind_onshore', 'solar_pv', # Renewables
'large_hydro', 'small_hydro', 'tidal_stream', 'shoreline_wave', 'Other Renewables',
'battery', 'Pumped Storage Hydroelectricity', 'pumped_hydro', # Storage
'H2', 'H2_turbine', # Hydrogen power
'CCGT', # Mid-merit gas
'Interconnectors Import', 'EU_import', # Imports
'coal', # Coal
'OCGT', # Peaker
'load_shedding' # Emergency (top)
]
# Group small contributors (<1% of total) into 'Other' categories
total_gen = p_by_carrier.sum().sum()
threshold = total_gen * 0.01 # 1% threshold
# Identify small bioenergy and renewable types to group
bioenergy_types = ['biogas', 'sewage_gas', 'advanced_biofuel', 'landfill_gas']
renewable_types = ['tidal_stream', 'shoreline_wave', 'small_hydro']
p_by_carrier_grouped = p_by_carrier.copy()
# Group small bioenergy
small_bio = [c for c in bioenergy_types if c in p_by_carrier.columns and p_by_carrier[c].sum() < threshold]
if len(small_bio) > 0:
p_by_carrier_grouped['Other Bioenergy'] = p_by_carrier_grouped[small_bio].sum(axis=1)
p_by_carrier_grouped = p_by_carrier_grouped.drop(columns=small_bio)
# Group small renewables
small_ren = [c for c in renewable_types if c in p_by_carrier.columns and p_by_carrier[c].sum() < threshold]
if len(small_ren) > 0:
p_by_carrier_grouped['Other Renewables'] = p_by_carrier_grouped[small_ren].sum(axis=1)
p_by_carrier_grouped = p_by_carrier_grouped.drop(columns=small_ren)
# Select columns with generation > 0
cols_with_gen = p_by_carrier_grouped.sum()[p_by_carrier_grouped.sum() > 0].index.tolist()
# Sort by stacking order
cols_to_plot = [c for c in stacking_order if c in cols_with_gen]
# Add any remaining columns not in stacking_order
cols_to_plot += [c for c in cols_with_gen if c not in cols_to_plot]
# Create interactive stacked area chart
fig = go.Figure()
for col in cols_to_plot: # Stack in order
fig.add_trace(go.Scatter(
x=p_by_carrier_grouped.index,
y=p_by_carrier_grouped[col] / 1e3, # Convert to GW
name=col.replace('_', ' ').title(),
mode='lines',
stackgroup='one',
fillcolor=color_map.get(col, '#CCCCCC'),
line=dict(width=0.5, color=color_map.get(col, '#CCCCCC')),
hovertemplate='<b>%{fullData.name}</b><br>' +
'Time: %{x}<br>' +
'Power: %{y:.2f} GW<br>' +
'<extra></extra>'
))
fig.update_layout(
title=dict(text='Generation Mix - Historical_2023_etys', x=0.5, xanchor='center'),
xaxis_title='Time',
yaxis_title='Generation (GW)',
hovermode='x unified',
height=600,
legend=dict(
orientation='v',
yanchor='top',
y=1,
xanchor='left',
x=1.02
),
template='plotly_white'
)
fig.show()
print('Stacking order: Baseload (bottom) -> Renewables -> Storage -> Gas -> Peakers (top)')
print('Tip: Click legend items to show/hide technologies. Double-click to isolate one.')
Stacking order: Baseload (bottom) -> Renewables -> Storage -> Gas -> Peakers (top)
Tip: Click legend items to show/hide technologies. Double-click to isolate one.
Interactive Network Topology#
Explore the network topology with interactive maps showing bus locations, transmission lines/links, and generators overlaid on a GB map.
Understanding Network Topology#
The maps below show the physical layout of Great Britain’s electricity infrastructure:
Buses: Connection points (substations, generation sites) shown as orange markers
Lines/Links: Transmission corridors carrying power between buses
Geography: Real GB coordinates using British National Grid projection
The network model captures how electricity flows are constrained by transmission capacity, creating locational price differences and potential congestion.
[5]:
# Interactive network map using PyPSA's explore() function
# Select a representative snapshot for visualization
snapshot_idx = len(network.snapshots) // 2
snapshot = network.snapshots[snapshot_idx]
print(f'Visualizing network at snapshot: {snapshot}')
# Detect network type based on bus count
is_zonal = len(network.buses) < 50 and len(network.lines) == 0 and len(network.links) > 0
is_reduced = len(network.buses) < 100 and len(network.buses) >= 30
network_type = 'Zonal' if is_zonal else ('Reduced' if is_reduced else 'ETYS')
print(f'Network type detected: {network_type} ({len(network.buses)} buses)')
# Network topology visualization (electricity only, no hydrogen)
try:
import copy
import folium
if 'x' in getattr(network.buses, 'columns', []) and 'y' in getattr(network.buses, 'columns', []):
# Filter out hydrogen buses and links
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas', 'H2 pipeline']
h2_buses = set()
if 'carrier' in network.buses.columns:
h2_buses = set(network.buses[network.buses.carrier.isin(h2_carriers)].index)
# Create filtered network copy
network_plot = copy.deepcopy(network)
# Remove hydrogen buses if any found
if len(h2_buses) > 0:
for bus_id in h2_buses:
if bus_id in network_plot.buses.index:
network_plot.remove('Bus', bus_id)
# Also filter hydrogen links
if len(network_plot.links) > 0 and 'carrier' in network_plot.links.columns:
h2_links = network_plot.links[network_plot.links.carrier.isin(h2_carriers)].index
for link_id in h2_links:
if link_id in network_plot.links.index:
network_plot.remove('Link', link_id)
# Convert OSGB -> WGS84 lon/lat for mapping
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
lon, lat = t.transform(network_plot.buses['x'].to_numpy(), network_plot.buses['y'].to_numpy())
network_plot.buses['x'] = lon
network_plot.buses['y'] = lat
print(f'Electricity network: {len(network_plot.buses)} buses, {len(network_plot.lines)} lines, {len(network_plot.links)} links')
# Create interactive map
m = network_plot.plot.explore(map_style='light', tooltip=True)
display(m)
print('✓ Network topology displayed (electricity only)')
else:
print('⚠️ No bus coordinates (x, y) found. Network visualization requires coordinate data.')
except Exception as e:
print(f'⚠️ Network topology map unavailable: {e}')
import traceback
traceback.print_exc()
Visualizing network at snapshot: 2023-05-18 12:00:00
Network type detected: ETYS (2044 buses)
Electricity network: 2044 buses, 1592 lines, 13 links
✓ Network topology displayed (electricity only)
[6]:
# Interactive network map showing transmission loading using PyPSA explore
try:
import copy
import folium
from folium import plugins
# Detect transmission element type (electricity only)
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = False
if has_lines:
# ETYS/Reduced network - use lines
loading = (network.lines_t.p0.loc[snapshot].abs() / network.lines.s_nom).fillna(0)
component_type = 'lines'
elif len(network.links) > 0:
# Zonal network - use internal electricity links only
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas', 'H2 pipeline']
is_external = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
is_h2_link = network.links.carrier.isin(h2_carriers)
internal_links = network.links[~is_external & ~is_h2_link]
if len(internal_links) > 0 and snapshot in network.links_t.p0.index:
link_flows = network.links_t.p0.loc[snapshot, internal_links.index].abs()
loading = (link_flows / internal_links.p_nom).fillna(0)
component_type = 'links'
has_internal_links = True
else:
loading = pd.Series(dtype=float)
if len(loading) > 0:
# Convert bus coordinates
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)
# Create folium map
center_lat = bus_coords['lat'].mean()
center_lon = bus_coords['lon'].mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')
# Color function
def get_color_by_loading(load_val):
if load_val >= 0.95: return '#DC143C' # Crimson (>95%)
elif load_val >= 0.80: return '#FF8C00' # Orange (80-95%)
elif load_val >= 0.60: return '#FFD700' # Gold (60-80%)
else: return '#228B22' # Green (<60%)
# Draw transmission elements with loading colors
if has_lines:
for line_id, load_val in loading.items():
line = network.lines.loc[line_id]
bus0_coords = [bus_coords.loc[line.bus0, 'lat'], bus_coords.loc[line.bus0, 'lon']]
bus1_coords = [bus_coords.loc[line.bus1, 'lat'], bus_coords.loc[line.bus1, 'lon']]
color = get_color_by_loading(load_val)
weight = 0.5 + (load_val * 4) # Width scales with loading
folium.PolyLine(
[bus0_coords, bus1_coords],
color=color,
weight=weight,
opacity=0.7,
tooltip=f'{line_id}: {load_val*100:.1f}% loaded'
).add_to(m)
elif has_internal_links:
for link_id, load_val in loading.items():
link = network.links.loc[link_id]
bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
color = get_color_by_loading(load_val)
weight = 0.5 + (load_val * 4)
folium.PolyLine(
[bus0_coords, bus1_coords],
color=color,
weight=weight,
opacity=0.7,
tooltip=f'{link_id}: {load_val*100:.1f}% loaded'
).add_to(m)
# Add small bus markers
for bus_id, row in bus_coords.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=2,
color='orange',
fill=True,
fillOpacity=0.6,
tooltip=bus_id
).add_to(m)
display(m)
# Print statistics
print(f'\nTransmission Loading ({component_type}) at {snapshot}:')
print(f' Max: {loading.max()*100:.1f}% | Mean: {loading.mean()*100:.1f}%')
print(f' >95% loaded: {(loading >= 0.95).sum()} | >80% loaded: {(loading >= 0.80).sum()}')
print('Color legend: Green (<60%) → Gold (60-80%) → Orange (80-95%) → Red (>95%)')
else:
print(f'⚠️ No transmission loading data available at {snapshot}')
except Exception as e:
print(f'⚠️ Transmission loading map unavailable: {e}')
import traceback
traceback.print_exc()
Make this Notebook Trusted to load map: File -> Trust Notebook
Transmission Loading (lines) at 2023-05-18 12:00:00:
Max: 100.0% | Mean: 8.8%
>95% loaded: 8 | >80% loaded: 11
Color legend: Green (<60%) → Gold (60-80%) → Orange (80-95%) → Red (>95%)
[7]:
# Interactive network map showing generator dispatch (electricity only)
try:
import folium
# Get total generation per bus at snapshot (exclude hydrogen generators)
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2']
elec_gens = network.generators[~network.generators.carrier.isin(h2_carriers)]
gen_per_bus = network.generators_t.p.loc[snapshot, elec_gens.index].groupby(elec_gens.bus).sum()
gen_per_bus = gen_per_bus[gen_per_bus > 0] # Only buses with generation
if len(gen_per_bus) > 0 and 'x' in getattr(network.buses, 'columns', []):
# Convert bus coordinates
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)
# Create folium map
center_lat = bus_coords['lat'].mean()
center_lon = bus_coords['lon'].mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')
# Draw light gray transmission lines in background
if len(network.lines) > 0:
for line_id, line in network.lines.iterrows():
bus0_coords = [bus_coords.loc[line.bus0, 'lat'], bus_coords.loc[line.bus0, 'lon']]
bus1_coords = [bus_coords.loc[line.bus1, 'lat'], bus_coords.loc[line.bus1, 'lon']]
folium.PolyLine([bus0_coords, bus1_coords], color='lightgray', weight=0.5, opacity=0.3).add_to(m)
# Color and size encoding for generation
max_gen = gen_per_bus.max()
def color_by_generation(gen_val):
norm = gen_val / max_gen
if norm < 0.25: return '#FFE4E1' # Misty rose
elif norm < 0.50: return '#FF9999' # Light red
elif norm < 0.75: return '#FF6666' # Red
else: return '#CC0000' # Dark red
# Draw bus markers sized and colored by generation
for bus_id, gen_mw in gen_per_bus.items():
coords = bus_coords.loc[bus_id]
norm_gen = gen_mw / max_gen
radius = 3 + (norm_gen * 12) # Size 3-15 pixels
color = color_by_generation(gen_mw)
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=radius,
color=color,
fill=True,
fillColor=color,
fillOpacity=0.7,
tooltip=f'{bus_id}: {gen_mw:.1f} MW'
).add_to(m)
display(m)
# Statistics
print(f'\nGenerator Dispatch at {snapshot}:')
print(f' Total generation: {gen_per_bus.sum():,.0f} MW')
print(f' Active buses: {len(gen_per_bus[gen_per_bus > 0])}')
print(f' Max at single bus: {gen_per_bus.max():,.0f} MW')
print('Color legend: Light gray (no gen) → Light red → Dark red (peak generation)')
else:
print('⚠️ No generator dispatch data available at this snapshot')
except Exception as e:
print(f'⚠️ Generator dispatch map unavailable: {e}')
import traceback
traceback.print_exc()
Make this Notebook Trusted to load map: File -> Trust Notebook
Generator Dispatch at 2023-05-18 12:00:00:
Total generation: 33,422 MW
Active buses: 2040
Max at single bus: 3,978 MW
Color legend: Light gray (no gen) → Light red → Dark red (peak generation)
[8]:
# Interactive installed capacity by carrier
generators_p_nom = network.generators.p_nom.groupby(
network.generators.carrier).sum().sort_values(ascending=True)
# Remove small contributors
generators_p_nom = generators_p_nom[generators_p_nom > 50]
# Create interactive bar chart
fig = go.Figure()
fig.add_trace(go.Bar(
y=generators_p_nom.index,
x=generators_p_nom.values,
orientation='h',
marker=dict(
color=[color_map.get(c, '#CCCCCC') for c in generators_p_nom.index],
),
text=[f'{val:,.0f} MW' for val in generators_p_nom.values],
textposition='auto',
hovertemplate='<b>%{y}</b><br>Capacity: %{x:,.0f} MW<extra></extra>'
))
fig.update_layout(
title=dict(text='Generator Capacity by Technology - Historical_2023_etys', x=0.5, xanchor='center'),
xaxis_title='Installed Capacity (MW)',
yaxis_title='',
height=500,
template='plotly_white',
showlegend=False
)
fig.show()
Storage Analysis#
[9]:
# Interactive storage dispatch and state of charge visualization
if len(network.storage_units) > 0:
p_storage = network.storage_units_t.p.sum(axis=1)
state_of_charge = network.storage_units_t.state_of_charge.sum(axis=1)
# Create subplot with two y-axes
fig = make_subplots(
rows=2, cols=1,
subplot_titles=('Storage Dispatch', 'State of Charge'),
vertical_spacing=0.12
)
# Storage dispatch
fig.add_trace(
go.Scatter(
x=p_storage.index,
y=p_storage.values,
name='Storage Dispatch',
line=dict(color='#32CD32', width=1.5),
fill='tozeroy',
fillcolor='rgba(50, 205, 50, 0.3)',
hovertemplate='Time: %{x}<br>Power: %{y:.0f} MW<extra></extra>'
),
row=1, col=1
)
# State of charge
fig.add_trace(
go.Scatter(
x=state_of_charge.index,
y=state_of_charge.values,
name='State of Charge',
line=dict(color='#00CED1', width=1.5),
fill='tozeroy',
fillcolor='rgba(0, 206, 209, 0.3)',
hovertemplate='Time: %{x}<br>Energy: %{y:.0f} MWh<extra></extra>'
),
row=2, col=1
)
fig.update_xaxes(title_text='Time', row=2, col=1)
fig.update_yaxes(title_text='Power (MW)', row=1, col=1)
fig.update_yaxes(title_text='Energy (MWh)', row=2, col=1)
fig.update_layout(
title=dict(text='Storage Analysis - Historical_2023_etys', x=0.5, xanchor='center'),
height=700,
hovermode='x unified',
template='plotly_white',
showlegend=True
)
fig.show()
# Storage statistics
print(f'\nStorage Statistics:')
print(f' Total storage capacity: {network.storage_units.p_nom.sum():,.0f} MW')
print(f' Total energy capacity: {network.storage_units.max_hours.sum() * network.storage_units.p_nom.sum():,.0f} MWh')
print(f' Peak discharge: {p_storage.max():,.0f} MW')
print(f' Peak charge: {p_storage.min():,.0f} MW')
print(f' Max state of charge: {state_of_charge.max():,.0f} MWh')
else:
print('No storage units in network')
Storage Statistics:
Total storage capacity: 5,216 MW
Total energy capacity: 970,234 MWh
Peak discharge: 4,366 MW
Peak charge: -3,978 MW
Max state of charge: 24,346 MWh
Hydrogen System Analysis#
The hydrogen system models the Power-to-Gas-to-Power pathway:
Electrolysis: Converts electricity to hydrogen (stored in central H2 bus)
H2 Storage: Central hydrogen storage (copper-plate model)
H2 Turbines: Converts hydrogen back to electricity when needed
About the Hydrogen System#
The hydrogen subsystem models the Power-to-Gas-to-Power pathway:
Electrolysis: Converts surplus renewable electricity to hydrogen (dashed turquoise links)
H2 Storage: Stores hydrogen for later use (purple markers)
H2 Turbines: Converts hydrogen back to electricity during high demand (red markers)
This provides seasonal energy storage, shifting summer renewable surplus to winter demand peaks.
[10]:
# Hydrogen System Analysis
# Check if hydrogen system exists
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
has_h2_stores = len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values
if has_electrolysis or has_h2_turbines:
print('=' * 80)
print('HYDROGEN SYSTEM SUMMARY')
print('=' * 80)
# Electrolysis
if has_electrolysis:
electrolysis = network.links[network.links.carrier == 'electrolysis']
elec_capacity = electrolysis.p_nom.sum()
elec_consumed = network.links_t.p0[electrolysis.index].sum().sum() / 1000 # GWh
elec_efficiency = electrolysis.efficiency.mean()
h2_produced = elec_consumed * elec_efficiency # GWh H2
print(f'\nELECTROLYSIS:')
print(f' Capacity: {elec_capacity:,.0f} MW ({len(electrolysis)} units)')
print(f' Efficiency: {elec_efficiency*100:.0f}%')
print(f' Electricity consumed: {elec_consumed:,.1f} GWh')
print(f' Hydrogen produced: {h2_produced:,.1f} GWh_H2')
# H2 Turbines
if has_h2_turbines:
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
turbine_capacity = h2_turbines.p_nom.sum()
h2_consumed = network.links_t.p0[h2_turbines.index].sum().sum() / 1000 # GWh H2
turbine_efficiency = h2_turbines.efficiency.mean()
elec_generated = h2_consumed * turbine_efficiency # GWh electricity
print(f'\nH2 TURBINES:')
print(f' Capacity: {turbine_capacity:,.0f} MW ({len(h2_turbines)} units)')
print(f' Efficiency: {turbine_efficiency*100:.0f}%')
print(f' Hydrogen consumed: {h2_consumed:,.1f} GWh_H2')
print(f' Electricity generated: {elec_generated:,.1f} GWh')
# H2 Storage
if has_h2_stores:
h2_stores = network.stores[network.stores.carrier == 'H2_gas']
h2_storage_capacity = h2_stores.e_nom.sum()
print(f'\nH2 STORAGE:')
print(f' Energy capacity: {h2_storage_capacity:,.0f} MWh ({h2_storage_capacity/1000:.0f} GWh)')
if h2_stores.index[0] in network.stores_t.e.columns:
soc = network.stores_t.e[h2_stores.index[0]]
print(f' State of charge range: {soc.min():,.0f} - {soc.max():,.0f} MWh')
print(f' Final state of charge: {soc.iloc[-1]:,.0f} MWh')
# Round-trip efficiency
if has_electrolysis and has_h2_turbines:
roundtrip_eff = elec_efficiency * turbine_efficiency
print(f'\nSYSTEM METRICS:')
print(f' Round-trip efficiency: {roundtrip_eff*100:.1f}%')
if elec_consumed > 0:
actual_roundtrip = elec_generated / elec_consumed
print(f' Actual energy ratio (elec out / elec in): {actual_roundtrip*100:.1f}%')
else:
print('No hydrogen system components found in this network.')
print('(Hydrogen components are added for future scenarios with H2 power generation)')
No hydrogen system components found in this network.
(Hydrogen components are added for future scenarios with H2 power generation)
[11]:
# Hydrogen System Dispatch Visualization
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
has_h2_stores = len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values
if has_electrolysis and has_h2_turbines:
electrolysis = network.links[network.links.carrier == 'electrolysis']
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
# Electrolysis power consumption (positive = consuming electricity)
elec_power = network.links_t.p0[electrolysis.index].sum(axis=1) / 1000 # GW
# H2 turbine power generation (p1 is output, negative means generation)
turbine_power = -network.links_t.p1[h2_turbines.index].sum(axis=1) / 1000 # GW
# Net hydrogen flow (positive = production, negative = consumption)
h2_production = elec_power * electrolysis.efficiency.mean() # GW H2 produced
h2_consumption = network.links_t.p0[h2_turbines.index].sum(axis=1) / 1000 # GW H2 consumed
# Create visualization
fig = make_subplots(
rows=3, cols=1,
subplot_titles=(
'Electrolysis (Electricity Consumed)',
'H2 Turbines (Electricity Generated)',
'H2 Storage State of Charge'
),
vertical_spacing=0.1
)
# Electrolysis power
fig.add_trace(
go.Scatter(
x=elec_power.index,
y=elec_power.values,
name='Electrolysis',
line=dict(color='#00CED1', width=1.5),
fill='tozeroy',
fillcolor='rgba(0, 206, 209, 0.3)',
hovertemplate='Time: %{x}<br>Power: %{y:.2f} GW<extra></extra>'
),
row=1, col=1
)
# H2 turbine power
fig.add_trace(
go.Scatter(
x=turbine_power.index,
y=turbine_power.values,
name='H2 Turbines',
line=dict(color='#FF6347', width=1.5),
fill='tozeroy',
fillcolor='rgba(255, 99, 71, 0.3)',
hovertemplate='Time: %{x}<br>Power: %{y:.2f} GW<extra></extra>'
),
row=2, col=1
)
# H2 storage state of charge
if has_h2_stores:
h2_stores = network.stores[network.stores.carrier == 'H2_gas']
if h2_stores.index[0] in network.stores_t.e.columns:
h2_soc = network.stores_t.e[h2_stores.index[0]] / 1000 # GWh
fig.add_trace(
go.Scatter(
x=h2_soc.index,
y=h2_soc.values,
name='H2 Storage',
line=dict(color='#9467BD', width=1.5),
fill='tozeroy',
fillcolor='rgba(148, 103, 189, 0.3)',
hovertemplate='Time: %{x}<br>Energy: %{y:.1f} GWh<extra></extra>'
),
row=3, col=1
)
fig.update_yaxes(title_text='Power (GW)', row=1, col=1)
fig.update_yaxes(title_text='Power (GW)', row=2, col=1)
fig.update_yaxes(title_text='Energy (GWh)', row=3, col=1)
fig.update_xaxes(title_text='Time', row=3, col=1)
fig.update_layout(
title=dict(text='Hydrogen System Dispatch - Historical_2023_etys', x=0.5, xanchor='center'),
height=900,
hovermode='x unified',
template='plotly_white',
showlegend=True
)
fig.show()
# Correlation analysis
print('\nDispatch Patterns:')
print(f' Peak electrolysis: {elec_power.max():.2f} GW at {elec_power.idxmax()}')
print(f' Peak H2 turbine: {turbine_power.max():.2f} GW at {turbine_power.idxmax()}')
print(f' Capacity factor (electrolysis): {elec_power.mean() / (electrolysis.p_nom.sum()/1000) * 100:.1f}%')
print(f' Capacity factor (H2 turbines): {turbine_power.mean() / (h2_turbines.p_nom.sum()/1000) * 100:.1f}%')
else:
print('Hydrogen dispatch visualization not available (no hydrogen system components)')
Hydrogen dispatch visualization not available (no hydrogen system components)
[12]:
# Interactive spatial map of hydrogen infrastructure
try:
import folium
has_h2_system = (
(len(network.links) > 0 and any(network.links.carrier.isin(['electrolysis', 'H2_turbine', 'H2 pipeline']))) or
(len(network.stores) > 0 and 'H2_gas' in network.stores.carrier.values)
)
if has_h2_system and 'x' in network.buses.columns:
# Convert all bus coordinates
t = Transformer.from_crs('EPSG:27700', 'EPSG:4326', always_xy=True)
bus_lon, bus_lat = t.transform(network.buses['x'].to_numpy(), network.buses['y'].to_numpy())
bus_coords = pd.DataFrame({'lon': bus_lon, 'lat': bus_lat}, index=network.buses.index)
# Create map centered on GB
center_lat = bus_coords['lat'].mean()
center_lon = bus_coords['lon'].mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=5, tiles='CartoDB positron')
# 1. Draw H2 pipelines (H2 pipeline links)
if len(network.links) > 0:
h2_pipelines = network.links[network.links.carrier == 'H2 pipeline']
if len(h2_pipelines) > 0:
for link_id, link in h2_pipelines.iterrows():
if link.bus0 in bus_coords.index and link.bus1 in bus_coords.index:
bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
capacity_mw = link.p_nom
folium.PolyLine(
[bus0_coords, bus1_coords],
color='#00CED1', # Cyan for H2 pipelines
weight=3,
opacity=0.7,
tooltip=f'H2 Pipeline: {link_id}<br>Capacity: {capacity_mw:.0f} MW'
).add_to(m)
print(f'Added {len(h2_pipelines)} H2 pipelines to map')
# 2. Draw electrolysis links (electricity -> H2) and facilities
if len(network.links) > 0:
electrolysis = network.links[network.links.carrier == 'electrolysis']
if len(electrolysis) > 0:
# Draw electrolysis links connecting electricity bus to H2 storage
for link_id, link in electrolysis.iterrows():
# bus0 is electricity input, bus1 is H2 output/storage
if link.bus0 in bus_coords.index and link.bus1 in bus_coords.index:
bus0_coords = [bus_coords.loc[link.bus0, 'lat'], bus_coords.loc[link.bus0, 'lon']]
bus1_coords = [bus_coords.loc[link.bus1, 'lat'], bus_coords.loc[link.bus1, 'lon']]
folium.PolyLine(
[bus0_coords, bus1_coords],
color='#40E0D0', # Turquoise for electrolysis links
weight=2,
opacity=0.6,
dash_array='5, 5', # Dashed line
tooltip=f'Electrolysis Link: {link_id}<br>Capacity: {link.p_nom:.0f} MW'
).add_to(m)
# Draw electrolysis facility markers at electricity input bus
for link_id, link in electrolysis.iterrows():
if link.bus0 in bus_coords.index:
coords = bus_coords.loc[link.bus0]
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=6 + (link.p_nom / 500), # Size by capacity
color='#40E0D0', # Turquoise
fill=True,
fillColor='#40E0D0',
fillOpacity=0.7,
tooltip=f'Electrolysis: {link_id}<br>Capacity: {link.p_nom:.0f} MW<br>Electricity bus: {link.bus0}<br>H2 bus: {link.bus1}'
).add_to(m)
print(f'Added {len(electrolysis)} electrolysis units and links to map')
# 3. Draw H2 turbine locations (H2 -> electricity)
if len(network.links) > 0:
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
if len(h2_turbines) > 0:
for link_id, link in h2_turbines.iterrows():
# bus1 is electricity output
if link.bus1 in bus_coords.index:
coords = bus_coords.loc[link.bus1]
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=6 + (link.p_nom / 500), # Size by capacity
color='#FF6347', # Tomato red
fill=True,
fillColor='#FF6347',
fillOpacity=0.7,
tooltip=f'H2 Turbine: {link_id}<br>Capacity: {link.p_nom:.0f} MW<br>Bus: {link.bus1}'
).add_to(m)
print(f'Added {len(h2_turbines)} H2 turbine units to map')
# 4. Draw H2 storage locations
if len(network.stores) > 0:
h2_stores = network.stores[network.stores.carrier == 'H2_gas']
if len(h2_stores) > 0:
for store_id, store in h2_stores.iterrows():
if store.bus in bus_coords.index:
coords = bus_coords.loc[store.bus]
energy_gwh = store.e_nom / 1000
folium.CircleMarker(
location=[coords['lat'], coords['lon']],
radius=8,
color='#9467BD', # Purple
fill=True,
fillColor='#9467BD',
fillOpacity=0.8,
tooltip=f'H2 Storage: {store_id}<br>Capacity: {energy_gwh:.1f} GWh<br>Bus: {store.bus}'
).add_to(m)
print(f'Added {len(h2_stores)} H2 storage facilities to map')
# Add legend
legend_html = '''<div style=\"position: fixed; bottom: 50px; right: 50px; width: 240px;
background-color: white; border:2px solid grey; z-index:9999; font-size:14px;
padding: 10px\">\n
<p style=\"margin:0; font-weight:bold;\">Hydrogen Infrastructure</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#00CED1; font-size:20px;\">━━</span> H2 Pipeline</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#40E0D0; font-size:12px;\">- - -</span> Electrolysis Link</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#40E0D0; font-size:20px;\">●</span> Electrolysis</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#FF6347; font-size:20px;\">●</span> H2 Turbine</p>\n
<p style=\"margin:5px 0;\"><span style=\"color:#9467BD; font-size:20px;\">●</span> H2 Storage</p>\n
</div>'''
m.get_root().html.add_child(folium.Element(legend_html))
display(m)
print('\\n✓ Hydrogen network map displayed')
print('Marker size indicates capacity (larger = higher capacity)')
else:
print('No hydrogen infrastructure found or missing bus coordinates')
except Exception as e:
print(f'⚠️ Hydrogen network map unavailable: {e}')
import traceback
traceback.print_exc()
No hydrogen infrastructure found or missing bus coordinates
Transmission Analysis#
[13]:
# Interactive transmission loading analysis with congestion hotspots
# Handles both lines (ETYS/Reduced) and links (Zonal)
snapshot_idx = len(network.snapshots) // 2
snapshot = network.snapshots[snapshot_idx]
# Determine which transmission elements to analyze
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = len(network.links) > 0
if has_lines:
# ETYS/Reduced network - use lines
loading = (network.lines_t.p0.loc[snapshot].abs() / network.lines.s_nom).fillna(0)
components = network.lines
element_type = 'Line'
elif has_internal_links:
# Zonal network - use internal links only (exclude external and hydrogen links)
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas']
is_external = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
is_h2_link = network.links.carrier.isin(h2_carriers)
internal_links = network.links[~is_external & ~is_h2_link]
if len(internal_links) > 0 and snapshot in network.links_t.p0.index:
link_flows = network.links_t.p0.loc[snapshot, internal_links.index].abs()
loading = (link_flows / internal_links.p_nom).fillna(0)
components = internal_links
element_type = 'Link'
else:
loading = pd.Series(dtype=float)
element_type = 'None'
else:
loading = pd.Series(dtype=float)
element_type = 'None'
if len(loading) > 0:
loading_sorted = loading.sort_values(ascending=False)
top_n = min(20, len(loading_sorted))
top_elements = loading_sorted.head(top_n)
# Create labels with bus names
element_labels = []
for elem_id in top_elements.index:
bus0 = components.loc[elem_id, 'bus0']
bus1 = components.loc[elem_id, 'bus1']
element_labels.append(f'{bus0} - {bus1}')
# Color code by congestion level
colors_list = [
'#DC143C' if x >= 0.95 else '#FF8C00' if x >= 0.8 else '#FFD700' if x >= 0.6 else '#228B22'
for x in top_elements.values
]
fig = go.Figure()
fig.add_trace(go.Bar(
y=element_labels,
x=top_elements.values * 100,
orientation='h',
marker=dict(color=colors_list),
text=[f'{val*100:.1f}%' for val in top_elements.values],
textposition='auto',
hovertemplate='<b>%{y}</b><br>Loading: %{x:.1f}%<extra></extra>'
))
# Add reference lines
fig.add_vline(x=95, line_dash='dash', line_color='red', annotation_text='95% (Congestion)')
fig.add_vline(x=100, line_dash='solid', line_color='darkred', annotation_text='100% (Limit)')
fig.update_layout(
title=dict(text=f'Top {top_n} Most Loaded {element_type}s - {scenario_id}', x=0.5, xanchor='center'),
xaxis_title=f'{element_type} Loading (%)',
yaxis_title='',
height=max(400, top_n * 25),
template='plotly_white',
showlegend=False
)
fig.show()
# Statistics
print(f'\n{element_type} Loading Statistics (snapshot: {snapshot})')
print(f' Total {element_type.lower()}s: {len(loading)}')
print(f' Max loading: {loading.max()*100:.1f}%')
print(f' Mean loading: {loading.mean()*100:.1f}%')
print(f' {element_type}s >95% loaded: {(loading >= 0.95).sum()}')
print(f' {element_type}s >80% loaded: {(loading >= 0.80).sum()}')
else:
print('No transmission elements found for loading analysis')
Line Loading Statistics (snapshot: 2023-05-18 12:00:00)
Total lines: 1592
Max loading: 100.0%
Mean loading: 8.8%
Lines >95% loaded: 8
Lines >80% loaded: 11
Renewable Curtailment Analysis#
[14]:
# Interactive renewable curtailment analysis
renewable_carriers = ['wind_onshore', 'wind_offshore', 'solar_pv']
# Create subplots for each renewable carrier
fig = make_subplots(
rows=len(renewable_carriers), cols=1,
subplot_titles=[carrier.replace('_', ' ').title() for carrier in renewable_carriers],
vertical_spacing=0.08
)
curtailment_stats = {}
for idx, carrier in enumerate(renewable_carriers, 1):
# Get generators of this carrier
gens = network.generators[network.generators.carrier == carrier]
if len(gens) == 0:
continue
# Calculate available and dispatched
p_nom_sum = gens.p_nom.sum()
# Check if any of these generators have p_max_pu time series data
# p_max_pu columns are indexed by generator name, not carrier
gens_with_pmax = [g for g in gens.index if g in network.generators_t.p_max_pu.columns]
if len(gens_with_pmax) > 0:
# Calculate available capacity from p_max_pu (weather-dependent availability)
available = network.generators_t.p_max_pu[gens_with_pmax].multiply(
gens.loc[gens_with_pmax, 'p_nom'], axis=1).sum(axis=1)
# Add any generators without p_max_pu at their full capacity
gens_without_pmax = [g for g in gens.index if g not in network.generators_t.p_max_pu.columns]
if len(gens_without_pmax) > 0:
available = available + gens.loc[gens_without_pmax, 'p_nom'].sum()
else:
# No time-varying availability - use installed capacity
available = pd.Series(p_nom_sum, index=network.snapshots)
dispatched = network.generators_t.p[gens.index].sum(axis=1)
curtailed = available - dispatched
# Add available capacity trace
fig.add_trace(
go.Scatter(
x=available.index,
y=available / 1000,
name='Available',
line=dict(color='orange', width=1.5),
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Available: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Add available capacity trace (orange line showing weather-dependent potential)
fig.add_trace(
go.Scatter(
x=available.index,
y=available / 1000,
name='Available',
line=dict(color='#FF8C00', width=1.5),
mode='lines',
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Available: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Add dispatched trace (filled area from zero)
fig.add_trace(
go.Scatter(
x=dispatched.index,
y=dispatched / 1000,
name='Dispatched',
line=dict(color='#228B22', width=1.5),
fill='tozeroy',
fillcolor='rgba(34, 139, 34, 0.6)',
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Dispatched: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Add curtailed as separate trace showing actual curtailment
curtailed_positive = curtailed.clip(lower=0) # Only show positive curtailment
fig.add_trace(
go.Scatter(
x=curtailed_positive.index,
y=curtailed_positive / 1000,
name='Curtailed',
line=dict(color='#DC143C', width=1.5, dash='dot'),
mode='lines',
legendgroup=carrier,
showlegend=(idx == 1),
hovertemplate='Curtailed: %{y:.2f} GW<extra></extra>'
),
row=idx, col=1
)
# Update y-axis label
fig.update_yaxes(title_text='Power (GW)', row=idx, col=1)
# Calculate statistics
curtailment_pct = (curtailed.sum() / available.sum() * 100) if available.sum() > 0 else 0
capacity_factor = (dispatched.sum() / (p_nom_sum * len(network.snapshots)) * 100)
curtailment_stats[carrier] = {
'capacity_mw': p_nom_sum,
'available_mwh': available.sum(),
'dispatched_mwh': dispatched.sum(),
'curtailed_mwh': curtailed.sum(),
'curtailment_pct': curtailment_pct,
'capacity_factor': capacity_factor
}
# Update layout
fig.update_xaxes(title_text='Time', row=len(renewable_carriers), col=1)
fig.update_layout(
title=dict(text='Renewable Curtailment Analysis - Historical_2023_etys', x=0.5, xanchor='center'),
height=900,
hovermode='x unified',
template='plotly_white'
)
fig.show()
# Print statistics
print('\nRenewable Curtailment Statistics:')
print('=' * 80)
for carrier, stats in curtailment_stats.items():
print(f"\n{carrier.replace('_', ' ').title()}:")
print(f" Installed Capacity: {stats['capacity_mw']:,.0f} MW")
print(f" Available Energy: {stats['available_mwh']:,.0f} MWh")
print(f" Dispatched Energy: {stats['dispatched_mwh']:,.0f} MWh")
print(f" Curtailed Energy: {stats['curtailed_mwh']:,.0f} MWh ({stats['curtailment_pct']:.1f}%)")
print(f" Capacity Factor: {stats['capacity_factor']:.1f}%")
Renewable Curtailment Statistics:
================================================================================
Wind Onshore:
Installed Capacity: 12,892 MW
Available Energy: 970,593 MWh
Dispatched Energy: 697,733 MWh
Curtailed Energy: 272,860 MWh (28.1%)
Capacity Factor: 32.2%
Wind Offshore:
Installed Capacity: 14,679 MW
Available Energy: 1,294,766 MWh
Dispatched Energy: 1,150,376 MWh
Curtailed Energy: 144,391 MWh (11.2%)
Capacity Factor: 46.6%
Solar Pv:
Installed Capacity: 9,019 MW
Available Energy: 280,556 MWh
Dispatched Energy: 276,071 MWh
Curtailed Energy: 4,485 MWh (1.6%)
Capacity Factor: 18.2%
Advanced Analysis#
Temporal Patterns and System Dynamics#
[15]:
# Interactive heatmap of hourly generation patterns
# Aggregate generation by hour of day and day of year
if len(network.snapshots) > 24:
total_gen = network.generators_t.p.sum(axis=1)
total_demand = network.loads_t.p_set.sum(axis=1)
# Create DataFrame with datetime index
df_temporal = pd.DataFrame({
'generation': total_gen,
'demand': total_demand,
'hour': total_gen.index.hour,
'day': total_gen.index.dayofyear
})
# Create hourly average pattern
hourly_pattern = df_temporal.groupby('hour')[['generation', 'demand']].mean()
# Interactive line plot of hourly patterns
fig = go.Figure()
fig.add_trace(go.Scatter(
x=hourly_pattern.index,
y=hourly_pattern['demand'] / 1000,
name='Average Demand',
line=dict(color='blue', width=2),
hovertemplate='Hour: %{x}<br>Demand: %{y:.1f} GW<extra></extra>'
))
fig.add_trace(go.Scatter(
x=hourly_pattern.index,
y=hourly_pattern['generation'] / 1000,
name='Average Generation',
line=dict(color='red', width=2),
hovertemplate='Hour: %{x}<br>Generation: %{y:.1f} GW<extra></extra>'
))
fig.update_layout(
title='Average Hourly Pattern',
xaxis_title='Hour of Day',
yaxis_title='Power (GW)',
template='plotly_white',
hovermode='x unified',
height=400
)
fig.show()
print(f'Peak demand hour: {hourly_pattern["demand"].idxmax()}:00')
print(f'Minimum demand hour: {hourly_pattern["demand"].idxmin()}:00')
else:
print('Insufficient timesteps for temporal pattern analysis')
Peak demand hour: 11:00
Minimum demand hour: 3:00
Congestion Hotspots Over Time#
Identify when and where transmission congestion occurs.
[16]:
# Calculate transmission loading over time (handles both lines and links)
has_lines = len(network.lines) > 0 and 's_nom' in network.lines.columns
has_internal_links = len(network.links) > 0
if has_lines:
# ETYS/Reduced network - use lines
loading_pct = (network.lines_t.p0.abs() / network.lines.s_nom * 100).fillna(0)
components = network.lines
element_type = 'Line'
elif has_internal_links:
# Zonal network - use internal links only (exclude external and hydrogen links)
h2_carriers = ['electrolysis', 'H2_turbine', 'H2_power', 'H2', 'H2_gas']
is_external = (
network.links.bus0.str.lower().str.contains('external', na=False) |
network.links.bus1.str.lower().str.contains('external', na=False)
)
is_h2_link = network.links.carrier.isin(h2_carriers)
internal_links = network.links[~is_external & ~is_h2_link]
if len(internal_links) > 0:
link_flows = network.links_t.p0[internal_links.index].abs()
loading_pct = (link_flows / internal_links.p_nom * 100).fillna(0)
components = internal_links
element_type = 'Link'
else:
loading_pct = pd.DataFrame()
element_type = 'None'
else:
loading_pct = pd.DataFrame()
element_type = 'None'
if len(loading_pct.columns) > 0:
# Count congested elements at each timestep (>95% loading)
congested_count = (loading_pct >= 95).sum(axis=1)
# Create time series plot
fig = go.Figure()
fig.add_trace(go.Scatter(
x=congested_count.index,
y=congested_count.values,
name=f'Congested {element_type}s',
line=dict(color='#DC143C', width=1.5),
fill='tozeroy',
fillcolor='rgba(220, 20, 60, 0.2)',
hovertemplate='Time: %{x}<br>Congested: %{y}<extra></extra>'
))
fig.update_layout(
title=f'Transmission Congestion Over Time ({element_type}s >95% Loaded)',
xaxis_title='Time',
yaxis_title=f'Number of Congested {element_type}s',
template='plotly_white',
height=400
)
fig.show()
# Find most frequently congested elements
congestion_frequency = (loading_pct >= 95).sum(axis=0).sort_values(ascending=False)
top_congested = congestion_frequency.head(10)
if len(top_congested) > 0 and top_congested.iloc[0] > 0:
print(f'\nMost Frequently Congested {element_type}s:')
print(f'{element_type}'.ljust(50) + ' Congested Timesteps')
print('=' * 75)
for elem_id, count in top_congested.items():
if count > 0:
bus0 = components.loc[elem_id, 'bus0']
bus1 = components.loc[elem_id, 'bus1']
pct = count / len(network.snapshots) * 100
print(f'{bus0} - {bus1:<40} {count:>4} ({pct:.1f}%)')
else:
print(f'\n✓ No significant transmission congestion detected')
else:
print('No transmission elements found for congestion analysis')
Most Frequently Congested Lines:
Line Congested Timesteps
===========================================================================
HAMB4A - STAH4A 117 (69.6%)
HARK13 - JUNV1- 106 (63.1%)
BOLN41 - NINF41 85 (50.6%)
WALP41 - RACO41 75 (44.6%)
FERO1- - FERO1S 71 (42.3%)
FERO1S - INVE1J 64 (38.1%)
BHLA11 - GLEN1Q 38 (22.6%)
GRNA1- - JUNV1- 14 (8.3%)
GLLU1Q - NETS1- 12 (7.1%)
INWI1Q - TORN1- 7 (4.2%)
Key Performance Indicators#
[17]:
# Calculate KPIs
print(f'SCENARIO: {scenario_id}')
print('=' * 80)
print()
total_demand = network.loads_t.p_set.sum().sum()
total_generation = network.generators_t.p.sum().sum()
# Load shedding
ls_gens = network.generators[network.generators.carrier == 'load_shedding']
load_shedding = network.generators_t.p[ls_gens.index].sum().sum() if len(ls_gens) > 0 else 0
print(f'DEMAND:')
print(f' Total demand: {total_demand:,.0f} MWh')
print(f' Load shedding: {load_shedding:,.0f} MWh ({load_shedding/total_demand*100:.2f}%)')
print(f' Demand satisfied: {total_demand - load_shedding:,.0f} MWh ({(total_demand - load_shedding)/total_demand*100:.2f}%)')
print()
# Renewable share
renewable_gens = network.generators[network.generators.carrier.isin(['wind_onshore', 'wind_offshore', 'solar_pv', 'large_hydro', 'small_hydro'])]
renewable_gen = network.generators_t.p[renewable_gens.index].sum().sum()
renewable_share = renewable_gen / (total_generation + load_shedding) * 100
print(f'GENERATION:')
print(f' Total generation: {total_generation:,.0f} MWh')
print(f' Renewable generation: {renewable_gen:,.0f} MWh')
print(f' Renewable share: {renewable_share:.1f}%')
print()
# System cost
if hasattr(network, 'objective'):
print(f'OPTIMIZATION:')
print(f' Total system cost: £{network.objective:,.2f}')
if total_generation + load_shedding > 0:
cost_per_mwh = network.objective / (total_generation + load_shedding)
print(f' Cost per MWh: £{cost_per_mwh:.2f}/MWh')
print()
# Line statistics
max_loading = abs(network.lines_t.p0 / network.lines.s_nom).max().max()
mean_loading = abs(network.lines_t.p0 / network.lines.s_nom).mean().mean()
n_congested = (abs(network.lines_t.p0 / network.lines.s_nom) >= 0.95).sum().sum()
print(f'TRANSMISSION:')
print(f' Total lines: {len(network.lines)}')
print(f' Max line loading: {max_loading*100:.1f}%')
print(f' Mean line loading: {mean_loading*100:.1f}%')
print(f' Timesteps with congestion (≥95%): {n_congested}')
print()
# Hydrogen system
has_electrolysis = len(network.links) > 0 and 'electrolysis' in network.links.carrier.values
has_h2_turbines = len(network.links) > 0 and 'H2_turbine' in network.links.carrier.values
if has_electrolysis or has_h2_turbines:
print(f'HYDROGEN SYSTEM:')
if has_electrolysis:
electrolysis = network.links[network.links.carrier == 'electrolysis']
elec_consumed = network.links_t.p0[electrolysis.index].sum().sum()
print(f' Electrolysis capacity: {electrolysis.p_nom.sum():,.0f} MW')
print(f' Electricity consumed: {elec_consumed:,.0f} MWh')
if has_h2_turbines:
h2_turbines = network.links[network.links.carrier == 'H2_turbine']
elec_generated = -network.links_t.p1[h2_turbines.index].sum().sum()
print(f' H2 turbine capacity: {h2_turbines.p_nom.sum():,.0f} MW')
print(f' Electricity generated: {elec_generated:,.0f} MWh')
if has_electrolysis and has_h2_turbines:
roundtrip = elec_generated / elec_consumed * 100 if elec_consumed > 0 else 0
print(f' Round-trip efficiency: {roundtrip:.1f}%')
print()
print('=' * 80)
SCENARIO: Historical_2023_etys
================================================================================
DEMAND:
Total demand: 4,805,482 MWh
Load shedding: 7 MWh (0.00%)
Demand satisfied: 4,805,475 MWh (100.00%)
GENERATION:
Total generation: 4,822,551 MWh
Renewable generation: 2,152,915 MWh
Renewable share: 44.6%
OPTIMIZATION:
Total system cost: £141,445,422.73
Cost per MWh: £29.33/MWh
TRANSMISSION:
Total lines: 1592
Max line loading: 100.0%
Mean line loading: 8.3%
Timesteps with congestion (≥95%): 631
================================================================================